function [AVGTimeVectorS]=MK008j(initN,ruleN,iters,steps,pa,pb,pc)
%
%copyright
%Ming Wu
%Dept.Computer Science and Engineering
%Michigan State University
%cite: PLOS ONE (2009) ˇ°A Dynamic Analysis of IRS-PKR Signaling in Liver Cells: A Discrete Modeling Approachˇ±
%
%Input: 
%initN:initial state vector, a state could be 0/1/2 or 9 (random)
%ruleN:state transition rule matrix, if i activates j, ruleN(i,j)=1, inhibition -1, otherwise 0
%iters:number of iterations, or number of cells in population
%steps:number of time-steps
%pa, pb, pc: to visualize time series of abitrary three network components
%
%Output: AVGTimeVectors: average states along the time steps for each component, each column for one component, each row for one time-step.
%

[tt,N]=size(initN);
CC=ones(1,N); CC=CC.*1.5;

AVGTimeVectorS=zeros(steps+1,N);
for it=1:iters
X0=initN; 

for i=1:N
    if X0(i)==9
%        X0(i)=(floor(3*rand()));
% different random initial distribution
        X0(i)=(round(2*rand()));
    end
end
%Adjustments/perturbations: silence/inhibition
%X0(9)=abs((round(2*rand())-1));
%X0(9)=2;
AVGTimeVector=zeros(steps+1,N);

AVGTimeVector(1,:)=X0;
endo=1;
for st=1:steps
    X1=ones(1,N); 
    %X1=-1*X1;
    rr=randperm(N);
    for j=1:N
        Act=0; Inhi=0;
        rrj=rr(j);
        for i=1:N
            if (ruleN(i,rrj)>0)&&(X0(i)==2)
                Act=Act+1;
            end
            if (ruleN(i,rrj)<0)&&(X0(i)==2)
                Inhi=Inhi+1;
            end
            if (ruleN(i,rrj)>0)&&(X0(i)==0)
                Act=Act-1;
            end
            if (ruleN(i,rrj)<0)&&(X0(i)==0)
                Inhi=Inhi-1;
            end                                  
        end
        if (rrj==5)&&(X1(4)<1)
            Act=0;inhi=0;
        end
            
        if Act>Inhi
            X1(rrj)=X1(rrj)+1;
            if X1(rrj)>2
                X1(rrj)=2;
            end
        else if Act<Inhi
                X1(rrj)=X1(rrj)-1;
                if X1(rrj)<0;
                    X1(rrj)=0;
                end
            end
        end
  %the first component is the stimulator    
       if rrj==1
           X1(rrj)=X0(1);
       end
       
 %Adjustments/perturbations: silence/inhibition 
 %      if rrj==3
 %          X1(rrj)=0;
 %      end
%       if rrj==4
%           X1(rrj)=0;
%       end
 %       if rrj==9
 %           X1(rrj)=X0(rrj);
 %       end

 %      if (rrj==1)&(st>30)
 %          X1(rrj)=0;
 %      end
 %      if (rrj==2)&(st>20)
 %          if X1(rrj)>1
 %              if (10*rand())<1.5;               
 %                  X1(rrj)=1;
 %              end
 %          else
 %              X1(rrj)=1;
 %          end
 %      end
%       if (rrj==2)&(st>25)
%               X1(rrj)=0.5;
%       end
%       if (rrj==2)&(st>40)
%                   X1(rrj)=0;
%       end               
                 
       X0(rrj)=X1(rrj);       
    end
    AVGTimeVector(st+1,:)=X1;
    X0=X1;
end
AVGTimeVectorS=AVGTimeVectorS+AVGTimeVector;
end
AVGTimeVectorS=AVGTimeVectorS./iters;

%to plot three time-series curve:
%XX=(1:steps+1);

%plot(XX,AVGTimeVectorS(:,pa),'g-',XX,AVGTimeVectorS(:,pb),'r-',XX,AVGTimeVectorS(:,pc),'b:')